.. _tutorial_2:
TUTORIAL 2 : NPT Lennard-Jones fluid
====================================
:Authors: Andrey V. Brukhno (andrey.brukhno{at}stfc.ac.uk), James Grant (r.j.grant{at}bath.ac.uk), and John Purton (john.purton{at}stfc.ac.uk)
Introduction
------------
In the NpT ensemble number of particles, pressure and temperature (*N,p,T*) are kept constant, which implies that apart from the particle moves the volume is also allowed to vary.
That is, an additional MC move, uniformly sampling the volume, is attempted with the acceptance probability:
.. math::
P_{\mathrm{acc}}([\mathbf{r}_{1},V_1] \rightarrow [\mathbf{r}_2,V_2]) = \min(1, \exp \{- \beta [U(\mathbf{r}_2) - U(\mathbf{r}_1) + p_{ext}(V_{2}-V_{1}) - N \beta^{-1} \ln(V_{2} / V_{1}) ] \} )
where :math:`p_{ext}` is the external pressure. The particle coordinates are assumed to be scaled accordingly.
NOTE: The last term under the exponent needs to be modified in the cases where, for example, linear dimensions of the (cubic) cell or :math:`\ln(V)` are directly sampled instead of volume (keywords **linear** and **log** respectively).
In order to proceed to NpT simulation, we will modify a copy of CONTROL file that we used in Tutorial 1 for NVT simulation, wherein we need to introduce the directives specifying the volume move and the accumulation of additional statistics.
Navigate to the subfolder **tutorial_2** (in the folder **exercises**). Our initial set-up contains duplicates of the input files for Tutorial 1. The CONFIG and FIELD files remain unchanged, so the initial configuration and interactions between particles will be identical to the NVT case. For your reference, the subdirectory **npt-dry** contains an example set of input and output files.
.. code-block:: html
:linenos:
NPT simulation of Lennard-Jones fluid
use ortho
finish
seeds 12 34 56 78 # Seed RNG seeds explicitly to the default
nbrlist auto # Use a neighbour list to speed up
# energy calculations
maxnonbondnbrs 512 # Maximum number of neighbours in neighbour list
temperature 1.4283461511745 # Corresponds to T*=1.1876;
# T(in K) = T* / BOLTZMAN
pressure 0.0179123655568
steps 110000 # Number of moves to perform in simulation
equilibration 10000 # Equilibration period: statistics
#are gathered after this period
print 10000 # Print statistics every 'print' moves
stack 10000 # Size of blocks for block averaging to obtain statistics
sample coord 10000 # How often to print configurations to ARCHIVE.000
yamldata 1000 # collect YAML stats every 1000 move
revconformat dlmonte # REVCON file is in DL_MONTE CONFIG format
archiveformat dlpoly4 # ARCHIVE.000/HISTORY.000/TRAJECTORY.000 format
# In this case: HISTORY.000 in DLPOLY4 style
move atom 1 512 # Move atoms with a weight of 512
LJ core
move volume cubic linear 1 # Move volume, box is cubic,
# linear scaling with a weight of 1
start
The directive to invoke volume moves is **move volume**. The parameters to specify in the present tutorial are **cubic**, **linear** and **1** for the weight of volume moves among other MC steps (see the manual for further details). Note that the weight of atom moves has been set to 512 (see **move atom** directive).
Volume moves are more computationally intensive than single atom moves, so the *rule of thumb* is to attempt one volume move each time every atom has been attempted to move (so called "pass" or "sweep" through the system).
You will need to add the lines::
sample volume 1 1.0 # sample volume every V-step, with the bin size of 1 A^3
move volume cubic linear 1 # Move volume, box is cubic,
# linear sampling with a weight of 1
and::
pressure 0.0179123655568
Exercise
--------
As before, in order to extract the volume sequence along the MC time from YAMLDATA.000 and plot it in gnuplot, use the commands::
[tutorial_2]$ strip_yaml.sh volume
[tutorial_2]$ gnuplot
gnuplot> plot './volume.dat' u 1:2 w l t "Volume(MC step)"
.. figure:: ./images/tutorial2-volume.jpg
:width: 640px
Questions to ask yourself:
^^^^^^^^^^^^^^^^^^^^^^^^^^
* Would it help equilibration in the *NpT* ensemle to start the simulation after *NVT* equilibration?
* What would be the effect of changing pressure and/or temperature?
* How to *restart* the simulation for continuation?
Try to increase the pressure and run a longer simulation starting with the last saved configuration (i.e. REVCON.000).
Further exercises
-----------------
Now try the extension exercises to learn more about functionality within DL_MONTE and to optimise your calculation:
:ref:`tut2_ex1`
.. :ref:`tut2_ex2`
.. Link to next tutorial
Or move on to :ref:`tutorial_3` and learn how to obtain radial distribution functions (RDF).